home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / mntlib44.zoo / mntlib / _addsubd.cpp next >
Encoding:
Text File  |  1993-06-17  |  7.2 KB  |  301 lines
  1. |
  2. | double floating point add/subtract routine
  3. |
  4. #ifndef __M68881__
  5.     .text
  6.     .even
  7.     .globl    __subdf3, ___subdf3
  8.     .globl    __adddf3, ___adddf3
  9. # ifndef    sfp004
  10. |
  11. | written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).
  12. | Based on a 80x86 floating point packet from comp.os.minix, written by P.Housel
  13. |
  14. | Revision 1.3.6 michal 05-93 (ntomczak@vm.ucs.ualberta.ca)
  15. |   + ensure that x - x always returns +0, says IEEE,
  16. |     unless x is Inf or NaN - then return NaN
  17. |
  18. | Revision 1.3.5 michal 05-93 (ntomczak@vm.ucs.ualberta.ca)
  19. |   + code smoothing
  20. |
  21. | patched by Olaf Flebbe (flebbe@tat.physik.uni-tuebingen.de)
  22. |
  23. | Revision 1.3.4 olaf 11-92 :
  24. |  + added support for NaN and infinities
  25. |    > floating point is now excellent!
  26. |
  27. |  -- still lacks trap handling for exceptions
  28. |  -- dont know the external representation of quiet and signaling NaN
  29. |     I decided 0x7fffffff,ffffffff to be a quiet NaN
  30. |     the rest should be signaling (but isnt)
  31. |
  32. | Revision 1.3.3 olaf 11-92 :
  33. |  + changed to get rid of rounding bits. a sticky register (d3) is
  34. |    sufficient.
  35. |
  36. | Revision 1.3.2 olaf 10-92 :
  37. |  + increased comparson by one again. (Dont understand, but it works)
  38. |  + corrected negation of rounding bits and mantissa
  39. |     >enquire now detects correct IEEE precision
  40. |     >paranoia now qualifies add/sub as correctly rounded
  41. |
  42. | Revision 1.3.1 olaf 10-92 :
  43. |  + increased comparison of exponents by one.
  44. |  + initialized sticky byte
  45. |  + corrected handling of rounding bits
  46. |     >paranoia now detects no SERIOUS DEFECTS any more
  47. |     ** Patches need _normdf Rev 1.6.1 (or higher) **
  48. |
  49. | Revision 1.3, kub 01-90 :
  50. | added support for denormalized numbers
  51. |
  52. | Revision 1.2, kub 01-90 :
  53. | replace far shifts by swaps to gain speed
  54. |
  55. | Revision 1.1, kub 12-89 :
  56. | Ported over to 68k assembler
  57. |
  58. | Revision 1.0:
  59. | original 8088 code from P.S.Housel
  60.  
  61. __subdf3:
  62. ___subdf3:
  63.     eorb    #0x80,sp@(12)    | reverse sign of v
  64. __adddf3:
  65. ___adddf3:
  66.     lea    sp@(4),a0    | pointer to u and v parameter
  67.     moveml    d2-d7,sp@-    | save registers
  68.     moveml    a0@,d4-d5/d6-d7    | d4-d5 = v, d6-d7 = u
  69.  
  70.     movel    #0x0fffff,d3
  71.     movel    d6,d0        | d0 = u.exp
  72.     andl    d3,d6        | remove exponent from u.mantissa
  73.     movel    d0,d2        | d2.h = u.sign
  74.     swap    d0
  75.     movew    d0,d2        | d2.l = u.sign
  76.  
  77.     movel    d4,d1        | d1 = v.exp
  78.     andl    d3,d4        | remove exponent from v.mantissa
  79.     swap    d1
  80.     eorw    d1,d2        | d2.l = u.sign ^ v.sign (in bit 15)
  81.     clrb    d2        | we will use the lowest byte as a flag
  82.     moveq    #15,d3
  83.     bclr    d3,d0        | kill sign bit u.exp
  84.     bclr    d3,d1        | kill sign bit v.exp
  85.     btst    d3,d2        | same sign for u and v?
  86.     beq    0f
  87.     cmpl    d0,d1        | different signs - maybe x - x ?
  88.     bne    0f
  89.     cmpl    d5,d7
  90.     seq    d2        | set 'cancellation' flag
  91.  
  92. 0:    lsrw    #4,d0        | keep here exponents only
  93.     lsrw    #4,d1
  94. |
  95. | Now perform testing of NaN and infinities
  96. |
  97.     movew    #0x7ff,d3
  98.     cmpw    d3,d0
  99.     beq    0f
  100.     cmpw    d3,d1
  101.     bne    nospec
  102.     bra    1f
  103. |
  104. |    u is special
  105. |
  106. 0:    tstb    d2        
  107.     bne    retnan        | cancellation of specials -> NaN
  108.     movel    d7,d0
  109.     orl    d6,d0
  110.     bne    retnan        | arith with Nan gives always NaN
  111.     addqw    #8,a0        | adding to an address propagates anyway
  112.     cmpw    d3,d1
  113.     bne    0f        | skip check for NaN if v not special
  114. |
  115. |    v is special
  116. |
  117. 1:    movel    d5,d0
  118.     orl    d4,d0
  119.     bne    retnan
  120. 0:    movel    a0@,d0        | copy infinity
  121.     moveq    #0,d1
  122.     bra     return
  123. |
  124. | return a quiet NaN
  125. |
  126. retnan: moveql    #-1,d1
  127.     movel    d1,d0
  128.     lsrl    #1,d0        | 0x7fffffff -> d0
  129.     bra    return
  130. |
  131. | Ok, no inifinty or NaN involved..
  132. |
  133. nospec: tstb    d2
  134.     beq    0f
  135.     moveq    #0,d0        | x - x hence we always return +0
  136.     movel    d0,d1
  137. return:    moveml    sp@+,d2-d7
  138.     rts
  139.  
  140. 0:    moveq    #20,d3
  141.     bset    d3,d6        | restore implied leading "1"
  142.     tstw    d0        | check for zero exponent - no leading "1"
  143.     bne    1f
  144.     bclr    d3,d6        | no implied leading "1", instead ...
  145.     addqw    #1,d0        | "normalize" exponent
  146. 1:
  147.     bset    d3,d4        | restore implied leading "1"
  148.     tstw    d1        | check for zero exponent - no leading "1"
  149.     bne    1f
  150.     bclr    d3,d4        | no implied leading "1", instead ...
  151.     addqw    #1,d1        | "normalize" exponent
  152. 1:
  153.     moveq    #0,d3        | init sticky register
  154.     negw    d1        | d1 = u.exp - v.exp
  155.     addw    d0,d1
  156.     beq    5f        | exponents are equal - no shifting neccessary
  157.     bgt    1f        | not equal but no exchange neccessary
  158.     exg    d4,d6        | exchange u and v
  159.     exg    d5,d7
  160.     subw    d1,d0        | d0 = u.exp - (u.exp - v.exp) = v.exp
  161.     negw    d1
  162.     tstw    d2        | d2.h = u.sign ^ (u.sign ^ v.sign) = v.sign
  163.     bpl    1f
  164.     bchg    #31,d2
  165. 1:
  166.     cmpw    #55,d1        | is u so much bigger that v is not
  167.     bge    7f        | significant ?
  168. |
  169. | shift mantissa left two digits, to allow cancellation of
  170. | most significant digit, while gaining an additional digit for
  171. | rounding.
  172. |
  173.     moveq    #1,d3
  174. 2:    addl    d7,d7
  175.     addxl    d6,d6
  176.     subqw    #1,d0        | decrement exponent
  177.     subqw    #1,d1        | decrement counter
  178.     dbeq    d3,2b
  179.     moveq    #0,d3
  180. |
  181. | now shift other mantissa right as fast as possible (almost).
  182. |
  183. 3:
  184.     cmpw    #16,d1        | see if fast rotate possible
  185.     blt    4f
  186.     orw    d5,d3        | set sticky word
  187.     movew    d4,d5        | rotate by swapping register halfs
  188.     swap    d5
  189.     clrw    d4
  190.     swap    d4
  191.     subqw    #8,d1
  192.     subqw    #8,d1
  193.     bra    3b
  194.  
  195. 0:    moveb   d5,d2        | use d2.b as scratch
  196.     andb    #1,d2        | test if 1 is shifted out
  197.     orb    d2,d3        | and put it in sticky
  198.     lsrl    #1,d4        | shift v.mant right the rest of the way
  199.     roxrl    #1,d5        | to line it up with u.mant
  200. 4:    dbra    d1,0b        | loop
  201.  
  202. 5:
  203.     tstw    d2        | are the signs equal ?
  204.     bpl    6f        | yes, no negate necessary
  205. |
  206. | negate second mantissa. One has to check the sticky word in order
  207. | to correct the twos complement.
  208. |
  209.     tstw    d3        |
  210.     beq     9f        | No correction necessary
  211.     moveq    #0,d1
  212.     addql   #1,d5
  213.     addxl   d1,d4
  214. 9:    negl    d5
  215.     negxl    d4
  216.  
  217. 6:
  218.     addl    d5,d7        | u.mant = u.mant + v.mant
  219.     addxl    d4,d6
  220.     bcs    7f        | need not negate
  221.     tstw    d2        | opposite signs ?
  222.     bpl    7f        | do not need to negate result
  223.  
  224.     negl    d7
  225.     negxl    d6
  226.     notl    d2        | switch sign
  227. 7:
  228.     movel    d6,d4        | move result for normalization
  229.     movel    d7,d5
  230.     moveq    #0,d1
  231.     tstl    d3
  232.     beq     8f
  233.     moveql   #-1,d1
  234. 8:    swap    d2        | put sign into d2 (exponent is in d0)
  235.     jmp    norm_df        | leave registers on stack for norm_df
  236. #else    sfp004
  237. | double precision floating point stuff for Atari-gcc using the SFP004
  238. | developed with gas
  239. |
  240. | double floating point add/subtract routine
  241. |
  242. | M. Ritzert (mjr at dmzrzu71)
  243. |
  244. | 4.10.1990
  245. |
  246. | no NAN checking implemented since the 68881 treats this situation "correct",
  247. | i.e. according to IEEE
  248.  
  249. | addresses of the 68881 data port. This choice is fastest when much data is
  250. | transferred between the two processors.
  251.  
  252. comm =     -6
  253. resp =    -16
  254. zahl =      0
  255.  
  256. | waiting loop ...
  257. |
  258. | wait:
  259. | ww:    cmpiw    #0x8900,a0@(resp)
  260. |     beq    ww
  261. | is coded directly by
  262. |    .long    0x0c688900, 0xfff067f8
  263.  
  264. __subdf3:
  265. ___subdf3:
  266. | double precision subtraction
  267. | sub second arg from fp0
  268.     lea    0xfffffa50:w,a0
  269.     movew    #0x5400,a0@(comm)    | load first argument to fp0
  270.     cmpiw    #0x8900,a0@(resp)    | check
  271.     movel    sp@(4),a0@
  272.     movel    sp@(8),a0@
  273.     movew    #0x5428,a0@(comm)
  274.     .long    0x0c688900, 0xfff067f8
  275.     movel    sp@(12),a0@
  276.     movel    sp@(16),a0@
  277.     movew    #0x7400,a0@(comm)    | result to d0/d1
  278.     .long    0x0c688900, 0xfff067f8
  279.     movel    a0@,d0
  280.     movel    a0@,d1
  281.      rts
  282.  
  283. __adddf3:
  284. ___adddf3:
  285.     lea    0xfffffa50:w,a0
  286.     movew    #0x5400,a0@(comm)        | load fp0
  287.     cmpiw    #0x8900,a0@(resp)        | got it?
  288.     movel    sp@(4),a0@            | take a hi from stack to FPU
  289.     movel    sp@(8),a0@            | take a lo from stack to FPU
  290.     movew    #0x5422,a0@(comm)        | add second arg to fp0
  291.     .long    0x0c688900, 0xfff067f8
  292.     movel    sp@(12),a0@            | move b hi from stack to FPU
  293.     movel    sp@(16),a0@            | move b lo from stack to FPU
  294.     movew    #0x7400,a0@(comm)        | result to d0/d1
  295.     .long    0x0c688900, 0xfff067f8
  296.     movel    a0@,d0                | download result
  297.     movel    a0@,d1                | download result
  298.      rts
  299. #endif    sfp004
  300. #endif    __M68881__
  301.